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ABSTRACT 

Although supernova explosions and stellar winds happen at very small scales, they affect the 
interstellar medium (ISM) at galactic scales and regulate the formation of a whole galaxy. Pre- 
vious attempts of mimicking these effects in simulations of galaxy formation use very simplified 
assumptions. We develop a much more realistic prescription for modeling the feedback, which 
minimizes any ad hoc sub-grid physics. We start with developing high resolution models of the 
ISM and formulate the conditions required for its realistic functionality: formation of multi-phase 
medium with hot chimneys, super-bubbles, cold molecular phase, and very slow consumption of 
gas. We find that this can be achieved only by doing what the real Universe does: formation 
of dense (> 10 H atoms cm" 3 ), cold (T rj 100 K) molecular phase, where the star formation 
happens, and which is disrupted by young stars. Another important ingredient is the runaway 
stars: massive binary stars ejected from molecular clouds when one of the companions becomes a 
supernova. Those stars can move to 10-100 parsecs away from molecular clouds before exploding 
themselves as supernovae. This greatly facilitates the feedback. Once those effects are imple- 
mented into cosmological simulations, galaxy formation proceeds more realistically. For example, 
we do not have the overcooling problem. The angular momentum problem (resulting in a too 
massive bulge) is also reduced substantially: the rotation curves are nearly flat. The galaxy 
formation also becomes more violent. Just as often observed in QSO absorption lines, there are 
substantial outflows from forming and active galaxies. At high redshifts we routinely find gas 
with few hundred km s _1 and occasionally 1000 — 2000 km s _1 . The gas has high metallicity, 
which may exceed the solar metallicity. The temperature of the gas in the outflows and in chim- 
neys can be very high: T — 10 7 — 10 8 K. The density profile of dark matter is still consistent 
with a cuspy profile. The simulations reproduce this picture only if the resolution is very high: 
better than 50 pc, which is 10 times better than the typical resolution in previous cosmological 
simulations. Our simulations of galaxy formation reach the resolution of 35 pc. 

Subject headings: hydrodynamics, methods: n-body simulations, ISM: general, stars: formation, galaxies: 
formation, evolution 



1. Introduction 

The current cosmological paradigm, the ACDM 
Universe, has successfully expl ained the overall as- 



semb l y of cosmic struc t ures (IB lumcnthal e t al . 
1984t bavis et al. 1 1 198.4 fSnergel et al. Il2007h . In 



this picture ordinary matter ( "baryons" ) , which 
emits and absorbs light, passively follows the evo- 
lution of the dark matter. This should be cor- 
rected, if we want to make a realistic theory of 
galaxy formation. It is necessary to include the 
physics of the gas and galaxy formation into the 



ACDM paradigm, because, after all, many obser- 
vational evidences of cosmic structures come from 
the light emitted by galaxies. 

Galaxy formation is driven by a complex set 
of physical processes with very different spatial 
scales. Radiative cooling, star formation and su- 
pernova explosions happen at scales less than 1 pc, 
but they affect the formation of a whole galaxy 
(|Dekel fc Silklll986f) . In addition, large-scale cos- 
mological processes, such as gas accretion through 
cosmic filaments, and galaxy mergers, control the 
galaxy assembly. As a result, a complex interplay 
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between very different processes drives the forma- 
tion of galaxies. Cosmological gasdynamical sim- 
ulations have become useful tools to study galaxy 
formation. 

In early cosmological simulations "galaxies" 

formed with too small disks and a s ignificant frac- 

tion of angular momentum was lost ([Navarro fe Steinmetz 
20001 ).^ The situation has improved in the last 



The problem is that other effects can be missed at 
the same time due to a lack of resolution and an 
inaccurate modeling of feedback. 

Another method is to introduce a sub-resolution 
model in which the energy from supernova ex- 
plosions is stored in an unresolved hot phase, 
which does not cool and loo ses energy through 



years. ISommcr-L arsen et al. I (120031 ) : lGovernato et al. 
(|2004) ; iKaufmann et al.l (|2006l) show that a sub- 



kpc resolution is necessary to prevent an arti- 
ficial loss of angular momentum. Recent im- 
provements in both the resolution and model- 
ing of the feedback have result ed in simulations 
with extended galactic disk s jGovernato et al 



2004; Robert son et al. 11200 4; Broo k et al 



Okamoto et al. 



Governato et al. 




200E 
12007 



D'Onghia et al. 



2004; 



2006; 



However, simulated 
galaxies are still too concentrated, and more real- 
istic simulations with better resolution and better 
physics are needed to reproduce the shape of the 
rotation curves of observed galaxies. 

Current simulations lack the necessary resolu- 
tion to follow correctly the effect of supernova ex- 
plosions in the ISM. Because of this lack of reso- 
lution, the modeling of stellar feedback has relied 
on ad-hoc assumptions about the effect of stel- 
lar feedback at scales unresolved by simulations 
(0.2-1 Kpc). Early attempts to introduce stellar 
feedback into simulations found the obstacle of a 
strong radiative cooling. The energy deposited by 
supernova explosions was quic kly radiated away 



without any effect in the ISM (jKatz Ill992h . Sev- 



eral shortcuts have been proposed to get around 
this over-cooling problem. 

The most common method is to artificially 
stop cooling when the st e llar energy is deposited 
( Gerritsen fc Icke 19971: Thacker fc Couchman 



2000 : Sommer-Larsen et al. I 2003 : iKeres et al. 
2005c iGovernato et al. 1120071 ) This approach pro 



longs the adiabatic phase of supernova explosion 
(the Sedov solution) to about 30 Myr. The moti- 
vation behind this ad-hoc assumption is that the 
combination of blastwaves from different super- 
nova explosions and turbulent motions produces 
hot bubbles much larger than individual super- 
nova remnants and last longer. All this effects are 
not resolved with the current resolution. They 
do not develop in a self-consistent way. Instead, 
the delay in the cooling is introduced by hand. 



the evaporation of cold cloud s (jYepes et al. 111997 



ISpringel fc Hernquist 20031 ). In this model, the 
only effect of stellar feedback is to regulate the star 
formation: the hot gas is coupled with the cold 
phase through cloud formation and evaporation. 
As a result, this high entropy gas is artificially 
trapped within the galactic disk. Thus, galactic 
winds are introduced in a simplified way in or- 
der to reproduce other natural effects of stellar 
feedback, such as galactic outflows. 

An alternative approach assu mes kinetic feed- 
back instead of thermal feedback ([Navarro fc White 
1993). In that case, the energy from supernova 



explosions or stellar winds is transfered to the ki- 
netic energy of the surrounding medium. This 
energy is not dissipated directly by radiative cool- 
ing. However, in order to resolve this effect accu- 
rately, simulations should be able to resolve the 
expansion of individual supernova explosions or 
the stellar winds from individual stars. Currently, 
this is not possible. At larger scales, the picture 
is more complicated. Different blastwaves from 
different supernova explosions can collide, dissi- 
pating their kinetic energy. The same dissipation 
of energy happens in collisions of stellar winds in 
stellar clusters. So, it is commonly assumed that 
most of the kinetic energy from stellar feedback 
is dissipated into thermal energy at the small- 
est scales resolved by simulations. Nevertheless, 
this feedback-heated gas can expand. As a result, 
thermal energy can be transfered to kinetic energy. 
The net results are flows at large scales powered 
by the thermal feedback. However, feedback heat- 
ing should dominates over radiative cooling: only 
in this case those flows are produced. 

To summarize, the main problems of current 
simulations of galaxy formation are the lack of the 
necessary resolution and too simplified models of 
the complex hydrodynamic processes in the mul- 
tiphase ISM. 

The galactic ISM has a very wide ra nge o f densi- 
ties and temperatur es (for review see Cox ( 20051) 
and iFerriere I ([200 ll )). Three distinct phases are 
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distinguished: the dense cold gas (giant molecu- 
lar clouds (GMC), cold HI gas or diffuse clouds) 
with densities above 10 cm -3 and temperatures 
bellow 100 K, the warm component with densi- 
ties between 0.1 and 1 cm~ 3 and temperatures 
of several thousands degrees, and the hot phase 
with temperatures above 10 5 K and densities bel- 
low 10~ 2 cm~ 3 . This multiphase medium is set 
by the competition of cooling and heating mech- 
anism and the onset of thermal instabilities. The 
hot ISM component (T > 10 5 K) is usually asso- 
ciated with gas heated by shocks. They can be 
produced by turbulent motions driven by gravi- 
tational and thermal instabilities. However, these 
turbule nt driven shocks can only heat the gas up to 
10 6 K (jWada fc Norman 1 120011 ). Only supernova 
explosions and st ellar winds can produce larger 
gas te mperatures (jMcCrav fe Snowlll979tlspitzer I 
1990h . 

2D and 3D hydrodynamical simulations of 
the ISM have enough resolution (parsecs) to re- 
solve the multi-phase nature of the ISM and 
to explore complicate d effects of stellar feed- 
back on different scales (|Rosen fe Bregman" 19951: 



back on dirterent scales ( Koscn &; hSreg] 
Scalo et al. 1 1 19981: iKorpi et al. II1999FT 



20001: Ide Avillez 



de Avillez 



Wada fe Norman 



\c Brcitschwerdt 
l200lL l2007t " lsivz et al 




There is much to learn from these simulations. 
However, they typically focus on conditions in the 
solar neighborhood, which are different from what 
one may expect during early stages of galaxy for- 
mation. Not always they follow the whole gas 
cycle: cooling, star formation, and stellar feed- 
back. For example, Ide Avillez fc Breitschwerdt 
(2004) include star formation but artificially re- 
strict the rate of supernova explosions around a 
fixed value. However, this rate could be much 
higher in large star forming regions. As a result, 
the effect of stellar feedback is underestimated in 
these regions. Nevertheless, the effect of the stel- 
lar feedback in the ISM, such as the formation of 
hot bubbles and super-bubbles is resolved. 

It is crucial to understand where and how 
the energy from massive stars is released back 
to the ISM. While a large fraction of massive 
stars are found in stellar clusters and OB as- 
sociations, 10-30% are found in the field, away 
from any molecula r cloud or stellar cluster (jGies I 



1987t IStone I 119911 ). This population have pecu- 



30 km s _1 , much higher than the velocity disper- 
sion of the po pulation of m assive stars in clusters 
(10 km s _1 ) (|Stone Ill99lh . Some of these stars 



have large peculiar veloc ities, up to 200 km s 1 
([Hoogerwerf et al~ l l2000h . This is why they are 
called runaway stars. The current scenario of the 
origin of runaway stars is the ejection of these 
massive stars from stellar clusters. There are two 
possible mechanisms of this ejection. One possi- 
bility is the ejection due t o a supernova explosion 



in a close binary system (jZwickv 1 11957c iBlaauw 
Il96ll ). The second mechanism is the ejection due 
to dynamical enco unters in the crowd ed regions 



of stellar clusters (jPoveda et al. 1119671 ). In spite 



of the fact that a significant fraction of the stellar 
feedback occurs far from star forming regions, no 
attention has been paid to its effect on the galaxy 
formation. 

We first study the effect of stellar feedback in 
the ISM, using simulations of a Kpc-scale piece of 
the ISM with few parsecs resolution. Then, we 
check if this picture holds when the resolution is 
degraded to the resolution that our cosmological 
simulations can achieve at high redshift. Finally, 
we study the effect of stellar feedback in galaxy 
formation at high redshift using cosmological hy- 
drodynamical simulations. 

This paper is organized as follows. Section 2 
describes the necessary conditions in which stellar 
feedback dominates over radiative cooling. Section 
3 describes all details of the modeling of stellar 
feedback. Section 4 shows Kpc-scale simulations 
of the ISM. Section 5 describes the cosmological 
simulations of galaxy formation. Finally, section 
6 is the discussion and conclusion. Throughout 
the paper we give quantities in physical units. 

2. Physical conditions for the heating 
regime 

The thermodynamical state of the gas depends 
on two competing processes: heating from stel- 
lar feedback and cooling from radiative processes. 
They appear as source and sink terms of internal 
energy in the equation of the first law of thermo- 
dynamics: 



du _ 



A 



(1) 



liar kinematics. Their velocity dispersion is about 



where u is the internal energy per unit volume, 
p is the pressure of the gas, and v is its velocity. 



3 



Parameter T is the heating rate due to stellar feed- 
back, and A is the net cooling rate from radiative 
processes. 

The heating rate from stellar feedback can be 
expressed as the rate of energy losses from a young 
and active single stellar population with a given 
density, p*, young : 



r' 



(2) 



where T' is the specific rate of energy losses of 
the stellar population according to its age. The 
cooling rate can be expressed as: 



A 



n%A>, 



(3) 



where nu is the hydrogen number density. 

2.1. Heating versus radiative cooling 

Now, we can ask ourselves under which con- 
ditions the feedback heating dominates over the 
radiative looses. Using the expression, nu = 
P ga s/(^Hm H ), where p gas is the gas density, n H 
is the molecular weight per hydrogen atom and 
m,H is the hydrogen mass, the condition for heat- 
ing (A < T) can be expressed as: 



n H A-' < 



P*, young 
Pgas 



(4) 



Using typical values, we can rewrite the condition 
for the heating regime in the following way: 



( HH ) 
V 0.1 cm- 3 / 



P*, young 



A' 



10 22 ergs 1 cm 3 



< 



(5) 



Pgas J Uo^ergs-iM- 1 



The cooling rate, A', is a strong function of gas 
temperature. So, the temperature and the density 
of the gas are two key properties in establishing 
the cooling or the heating regimes. The following 
two examples illustrate common situations. 

At temperatures around 10 4 K, the cooling rate 
is close to its maximum value. We use A' = 10~ 22 
erg s _1 cm 3 as a fiducial value. In this case 
eq.© shows that the heating overcomes the cool- 
ing only at very low densities njj < 0.1 cm -3 , op- 
timistically assuming that the ratio of densities, 
P*, young/ Pgas is about unity. As a result, stellar 



feedback is not able to heat the gas beyond 10 4 
K for densities higher than 0.1 cm~ 3 and typical 
values of 17". This is the well known overcooling 
problem for simulations, which allow cooling only 
to a temperature of 10 4 K at which the star for- 
mation is assumed to happen. The energy from 
stellar feedback is radiated away very efficiently 
and the thermal feedback cannot play any role. In 
this case one needs to invoke "subgrid physics" - 
a guess how the system should react to the energy 
released by the stars. 

The situation is completely different if the gas 
is allowed to cool to 100 K. The cooling is very 
inefficient at that temperature: A' = 10 -25 erg 
s -i cm 3 g Q ^ s t e ll ar feedback can produce the net 
gas heating even if the density is large: nn ~ 100 
cm -3 for p*, young ~ Pgas- Our conclusion is that 
simulations should include cooling process bellow 
10 4 K. The cold phase should be resolved in order 
to get a high efficiency of stellar feedback. 

However, heating to high temperatures is still 
problematic because as the gas is heated, the 
cooling rate increases. So, the peak of the cool- 
ing rate at 10 4 K is a bottle-neck for heating 
gas to higher temperatures. Nevertheless, tem- 
peratures of diffuse gas as high as 10 6 — 10 7 K 
have been observed around star -forming regions 
such as the Rosette n ebulae ( Townslev et al 



2003; IWang et al 
2003 ). and the Orion 
2005 



20071 ). M17 ( Townslev et al. 
nebula ( Feigelson et al. 



Guedel et al. 1 120071 ). The main question 



is how young and massive stars can heat their sur- 
rounding medium to these high temperatures, if 
the original medium, in which they were born had 
high densities. 

The answer likely depends on the distance from 
those young stellar clusters. At small 1-2 pc dis- 
tances it is likely to have the collisions of stel 



lar w inds ( Townslev et al. 20031 : Feigelson et al 



2005) . At larger distances the heating is related 



with the formation of superbubbles: the cumula- 
tive effect of winds and shocks generated by many 
young stars. One way or another, the density of 
gas around the young stellar population decreases 
and the ratio p*, VO ung/Pgas increases as the over- 
pressured bubble of gas expands. Once the density 
goes below 0.1 cm -3 , eq.© can be fulfilled even 
at 10 4 K. The net result is a heating regime, in 
which the surrounding gas can be heated to very 
high temperatures. In other words, the process 
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starts with the expanding bubbles at low temper- 
atures and then proceeds to a runaway overheating 
regime. 

As an example, we consider a typical GMC with 
a mass of 10 5 M Q and a size of 50 pc. These 
are the typical va lues found in recent ca talogs 
of GMCs in M33 (iRosolowskv et al J l2007h. M31 
and the Milky way (|Sheth et al.ll2008l ). Therefore, 
the mean density is njj — 50 cm~ 3 . This value 
seems low compared with typical observed den- 
sities of molecular clouds. However, GMCs are 
highly clumpy. High-density clumps are usually 
embedded in a low density inter-clump medium. 
As a result, the volume-averaged density inside 
clouds is much smaller than the typ ical observed 
mass- weighted density ( McKee 19991 ). 

Now, we consider an Orion-like stellar cluster 
formed at the center of the cloud. The mass of the 
stellar cluster is 5 x 10 3 M ( Hillenbrand k. Hartmanr] 
1998). In a region of mass 10 4 M , the stellar 
cluster has the ratio /0*, yO ung/Pgas equal to 0.5, 
and the condition for heating, eq.©, is fulfilled. 
This heating produces an over-pressured hot bub- 
ble with a pressure 100 times higher than the 
surrounding unperturbed medium. As a result, 
the bubble expands, the density decreases, and 
the ratio p*, y0 ung/Pgas increases. Then, we get a 
runawa y bubble, which pro ceeds to blowing away 
all gas (jKroupa et aDl200lh . 

Simulations should resolve the expansion of 
bubbles over-pressured by stellar feedback. The 
density of young (and active) stars and the den- 
sity of gas should be comparable at the smallest 
scales resolved by the simulations. The minimum 
value of the ratio p*, y0 ung/Pgas depends on the 
gas density (eq. [5]). For moderate gas densities, 
riff = 10 — 100 cm' 3 , the above ratio should be 
around 0.1-1. 

The above condition can be achieved if the star 
formation efficiency, the fraction of the progeni- 
tor cloud consumed in stars is 10%-50% at the 
resolution scale. This high efficiency is consis- 
tent with the observed val ue of 10%-40% found in 
Galactic stellar cluste r s, (Greene &; Yound 19921 : 
Elmegreen et al.ll200Ct Ikroupa et al.l l200lh . Due 
to the fact that 80% of th e Galactic star form ation 
occurs in stellar clusters ( Lada &: Ladal l2003). this 
high efficiency of star formation should be consid- 
ered in any star formation model which can resolve 
the sites where star formation occurs. 



2.2. Local gravity versus pressure gradient 

As we saw in the previous section, low densities 
are required in order to heat the gas beyond the 
peak of the cooling curve. Stellar feedback should 
evacuate the gas by creating an expanding bub- 
ble around young stellar clusters. However, the 
over-pressured bubble expands only if the pressure 
gradient overcomes self-gravity. 

If we consider an over-pressured bubble of ra- 
dius R in a homogeneous medium of density p, 
we can derive a Jeans-instability type of condi- 
tion. As a result, the bubble expands only if the 
difference in pressure with its surroundings, AP, 
satisfies the following relationship: 



4-7T 

AP/k > ^G(pR) 2 



10" 1 (n H Rp C 



(6) 



where k is the Boltzmann constant, G is the gravi- 
tational constant, and i? pc is the radius in pc. The 
above equation sets the conditions for the bubble 
expansion. For the Galact ic plane th e pressure 



is P/k ~ 2 x 10" cm- 3 if (|Cox 1 120051 ). For ex- 



ample, a region of 50 pc in radius and a density 
of 100 cm -3 will only expand, if the difference in 
pressure is bigger than 2 x 10 6 cm -3 K. This can be 
achieved, if the bubble is over-pressured by more 
than 100 times. Stellar feedback can produce this 
overpressure just by raising the temperature from 
100 K to 10 4 K. The resulted over-pressured re- 
gion will expand, and the density as well as the 
cooling rate will decrease. So, the efficiency of 
stellar feedback increases, raising the temperature 
and pressure further. 

Eq. [5] also sets a upper limit on the resolution. 
Using the equation of state of the ideal gas P = 
nkT, where n is the mean number density and T 
is the temperature of the gas, the over-pressured 
bubble should be resolved with a spatial resolution 
X pc = Rp C /2, such that the expansion is resolved: 



X. 



pc 



75pc) ~ ^10 4 k) VlOcm" 3 



( T 



-l 



(7) 



As a result, for typical values of these over- 
pressured regions, the resolution should be bet- 
ter than ~ 70 pc. Otherwise, the bubble cannot 
overcome its self-gravity and cannot expand. 
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3. Stellar feedback model 

We assume a model of thermal feedback for the 
injection of energy from stellar winds and super- 
nova explosions. The kinetic energy from these 
processes is efficiently dissipated into thermal en- 
ergy due to shocks at scales bellow the spatial res- 
olution. 

The net thermal rate (r — A) is used to update 
the internal energy in each step of the simulation. 
This approach is rather different than the deposi- 
tion of energy. Instead, the energy injection from 
stellar feedback is treated in a self-consistent way 
along with the radiative looses. 

3.1. Heating rate from stellar feedback 

The heating rate from stellar feedback in a 
given volume element is modeled as the rate of en- 
ergy losses from a set of single stellar populations 
present in that volume. This is just a generaliza- 
tion of eq.j2]): 



(8) 



where Mi and U are the mass and the age of each 
single stellar population. 

The modeling of the specific release of energy 
over time, T' , is motivated by the results from pop- 
ulation synthesis codes , such as STARBURST99 



( Leitherer et al. 1999). Figure [T] shows differ- 
ent models of V and the results of a STAR- 
BURST99 computation with a Miller-Scalo IMF 
from 0.1 M to 100 M . Parameter V is dom- 
inated by stellar winds from massive OB main- 
sequence stars and WR stars during the first 
few Myr. Later the energy is produced by core- 
collapse supernovae from stars more massive than 
8 Mq. After 40 Myr, the release of energy comes 
from stellar winds of AGB stars and other less 
powerful sources, and the injection rate drops 6 
orders of magnitudes. Supernovae la dominate the 
feedback at much longer time-scales. We assume 
a peak of the SNIa rate at 1 Gyr. However, this 
peak is 3 orders of magnitude lower than the con- 
tribution from core-collapse supernovae. This is 
because the energy from a population of SNIa is 
diluted over a much longer time scale than the en- 
ergy from core-collapse supernovae. 

We model V with a constant rate of 1.18 x 10 34 
erg s^M"" 1 over 40 Myr. This is equivalent to the 
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Fig. 1. — Rate of energy losses per unit mass from 
a single stellar population. Top panel shows the 
results from the STARBURST99 code, assuming a 
Miller-Scalo IMF for a mass range (0.1 - 100) Mq. 
The dotted line shows the contribution of super- 
nova explosions and the full line shows the total 
rate. Although supernova explosions dominate the 
overall energy release, stellar winds are the only 
mechanism of energy release during the first few 
Myr. Middle and bottom panels show two dif- 
ferent models: a constant feedback model and a 
model of stellar wind plus core-collapse supernova. 
Although the total energy released is the same 
in both models, the SN model is more elemen- 
tary and takes into account the explosive nature 
of core-collapse supernova. 

injection of 2 x 10 51 erg of energy from stellar winds 
and supernova explosions per each massive star 
with M > 8 M Q during its lifetime. We assume a 
Miller-Scalo IMF in the mass range (0.1-100) Mq. 
Note that this constant heating rate is the sum of 
the contributions from all massive stars in a single 
stellar population. We also consider a more sim- 
ple model, which we call a SN model. In this case 
10 51 ergs is injected at constant rate due to stellar 
winds over 10 or 40 Myr. Then it follows a strong 
peak of energy release due to the supernova explo- 
sion, in which 10 51 erg are released during 10 5 yrs 
- the typical age of young supernova remnants. 
Although the total energy released is the same in 
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both models, the SN model takes into account the 
explosive nature of core-collapse supernovae. 



metallicity-dependent cooling and UV heating due 
to a cosmological ionizing background. 



3.2. A model of runaway stars 



Results of ISM runs 



The effect of runaway stars is implemented by 
adding a random velocity to a fraction of stellar 
particles (10%-30%). This extra velocity has a 
random orientation and the value is taken from 
an exponential distribution with a characteristic 
scale of 17 km s - 1 . This choice is motiv ated 
by Hipparcos data (IHoogerwerf et al. 2000) and 



Monte-Carlo simulations (|Drav et al. 2005). For 



compar i son, a Gaussian distribution is also used 



([Stone 1 11991). However, the effect of runaway 



stars in the ISM is not very sensitive to the de- 
tails of this velocity distribution. 

3.3. Radiative cooling 

Radiative cooling counterbalances feedback 
heating. So it is very important to have an accu- 
rate model of radiative cooling in order to study 
the net effect of stellar feedback in the ISM. 

We use the mode l of radiative cooling described 
~~ ( 20031) . 



Kravtsov 



It is a metallicity-dependent 
cooling plus a UV hea ting due to a cosmolog i- 
cal ionizing background (jHaardt fc Madau 1 11996). 
The model includes Compton heating/cooling and 
molecular cooling. The temperature range of the 
model is between 10 2 K < T < 10 9 K. Thus, 
this model includes cooling below 10 4 K and the 
gas can reach the thermodynamical conditions of 
molecular clouds. As we saw in section 2, this is 
crucial for the efficiency of the stellar feedback. 

The cooling and heating rates from radiative 
processes are tabulated using the CLOUDY code 
(version 96b4; Ferland et al. 1998). As a result, 
the net cooling rate from radiative processes, A', is 
available for a given density, temperature, metal- 
licity and redshift. 

3.4. Description of the code 

The numerical simulations were performed us- 
ing the Eulerian gasdyna mics + N-body Adap- 



tivc Refinement Tree code (Kravtsov et al. 1997 



Kravtsov Il999l | 20031 ). The physical processes 
of the gas include star formation, stellar feed- 
back, metal enrichment, self-consistent advec- 
tion of metals, cooling and heating rates from 



Our first step in the understanding of stellar 
feedback in galaxies is to understand its effect in 
the ISM at galactic scales. Therefore, we run sim- 
ulations of a 4 x 4 x 4 Kpc 3 piece of a galactic 
disk with 8-16 pc resolution. These simulations 
fully resolve the effect of massive stars at galactic 
scales. So, resolution is not longer an issue. 

We can use this ISM-scale simulation as a 
benchmark for the effect of stellar feedback at 
galactic scales. Then, we can degrade the reso- 
lution to see which model of feedback reproduces 
the same overall picture at lower resolution. These 
simulations can then be used as testing grounds 
for these models at different resolutions. They tell 
us what are the necessary ingredients to reproduce 
the truly effect of stellar feedback at the resolution 
that we can afford in cosmological simulations of 
galaxy formation. 

We want to see the effect of stellar feedback in 
the typical conditions of normal disk galaxies with 
moderate gas surface densities. So, we are not 
modeling starburst galaxies with large amounts of 
gas and high star formation rates. This type of 
study will be done in the future. 

A 4 Kpc box of ISM represents a significant 
piece of a galactic disk. The simulation resolves 
the dense galactic plane, where molecular clouds 
are formed. This is important to follow star for- 
mation correctly. At the same time, the simula- 
tion follows the gas at few Kpc above the galactic 
plane. This height is similar to the scale-height of 
the diffuse warm phase of the ISM (Cox 2005). 

The simulation includes radiative cooling and 
UV heating from a uniform UV field at redshift 
as described in section 3. Star formation happens 
in the highest density peaks with a density thresh- 
old of 100 cm~ 3 . In each star formation event, 
5 % of the mass in gas inside a volume element 
is converted into a stellar particle with a mass 
of 88 M within a time-step set by the Courant 
condition (~ 2 xlO 3 yr). The supernova model 
was used for stellar feedback and SNIa was not 
included. The metallicity was assumed solar and 
constant throughout the simulation. 
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Fig. 2. — Formation of a galactic chimney. Edge- 
on slices through the simulation show density, 
temperature and velocity in the vertical direction, 
perpendicular to the galactic plane. The bottom 
panel shows gas column density. The chimney out- 
flow is not a homogeneous, coherent flow: it is tur- 
bulent and has dense and cold clumps embedded 
into the flow. The core of the chimney reaches 
10 7 -10 8 K. Outflow velocities exceed 10 3 km s _1 . 
This hot material is able to escape the disk and 
generate a galactic wind. 



4.1. Initial conditions 

The initial distribution of gas density is uniform 
in the x and y directions of the box. In the z- 
direction, the density profile declines at both sides 
of the middle plane, z = zq — 2 Kpc. This plane 
defines the galactic plane for this ISM model: 



no cosh 



Z - Zq 
Zd 



(9) 



where no is the gas density in that plane and Zd is 
the scale-height. 

The choice of parameters sets the conditions of 
a quiescent normal galactic disk, no = 1 cm -3 
and Zd — 250 pc. Thus, the surface density is 
, E gas = 16 M pc -2 . The system is originally 
in hydrostatic equilibrium with a temperature of 
10 4 K. No stars are present at the beginning of 
the simulation. The box has open boundaries in 
the z-direction. So, all material that cross these 
boundaries escapes the system. 

The initial velocity field consists of a sum of 
plane-parallel velocity waves: 



i x = 22 A x (i,j, k) sin(fc • r)exp- 

i,j,k 

■iy = A y (i,j, k) sin(fc • r) exp - 
u z = X! A z (i,j, k)sin(k ■ r)exp- 



Z~ Zq 
Zd 

Z - Zq 
Zd 

Z - Zq 
Zd 



(10) 



(11) 



(12) 



The amplitudes are taken from a Gaussian field 
with a tilted power spectrum, Pf, cx fc -3 , where k 
is the wavenumber, k — ^jj\/i 2 + j 2 + k 2 . i,j and 
k are integers running from -20 to 20 (excluding 
0) and ito = 20 km s _1 . This is a typic al spectrum 
of a compressible turbulent medium (jKraichnan 
Il967t IVazauez-Semadeni et al.lll995h . 



A x (i,j,k) = u Q — 



Rgo 



A y (i,j,k) = u Q — 
A z (i,j,k) = u 



(i 2 + f + fc 2 ) 3 / 2 

^Gauss 



(i 2 +f + k 2 ) 3 / 2 



(i 2 + f + fc 2 ) 3 / 2 



(13) 
(14) 
(15) 



Rgq 



is a random number taken from a Gaussian 



distribution. 
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Fig. 3. — Top panel: Star formation rate surface 
density of the whole simulation. The value shown 
is also averaged over a period of ~ 2 xlO 5 yr 
(100 time-steps). After an initial burst, the star 
form ation rate surface den sity is consistent with 
the iKennicutt et al. ( 20071 ) empirical fit (horizon- 
tal line) . Bottom panel: Fraction of volume filled 
with each gas phase over time. The volume oc- 
cupied by the warm and the hot phase oscillates. 
The hot phase dominates after a burst of star for- 
mation and the warm phase dominates when the 
gas is cooled down. The cold phase covers a small 
volume, which remains constant after an initial 
collapse. 



4.2. Galactic Chimney formation 

At the beginning of the simulation, the gas 
starts to move according to the turbulent veloc- 
ity field. As a result, the gas accumulates where 
different flows converge and molecular clouds 
naturally appear in form of filaments and shells. 
However, around 90% of the volume is filled with 
warm and diffuse gas heated by UV background. 
Star formation occurs in the cores of the cold 
phase. Newly formed massive stars inject energy 
and a cavity filled with hot and very diffuse gas is 
formed. This over-pressured material expands and 
the net result is the formation of super-bubbles. 



This hot gas cannot stay in the plane of the disk, as 
a result, the bubble expands faster in the direction 
perpendicular to the disk, because the density de- 
clines in that direc tion. The bubble develop s into 
a galactic chimney ( Norman fc Ikeuchll 19891 ) . The 
chimney outflow does not look as a homogeneous, 
coherent flow. Instead, the chimney is turbulent 
and has dense and cold clumps embedded into the 
flow. Eventually, the gas expands in the halo and 
cools (Figure^]). 

Another interesting feature seen in this model 
is a population of isolated bubbles in the warm 
medium. These are the results of individual su- 
pernova explosions of runaway stars 

4.3. Star formation rate 

After an initial burst of star formation, the 
star formation rate is nearly constant for the rest 
of the evolution (Figure [3]). We found a low 
star formation rate surface density, Ssfr = 3 x 
1O~ 3 M yr -1 Kpc~ 2 , temporally averaged over a 
period of 2 x 10 5 yr (100 time-steps). This value 
is consistent with the expected value from the 
correlation between the star formation rate sur- 
face density an d the gas surface density found in 
nearb y galaxies (jKennicutt 1998: IKennicutt et alj 
20071 ). For a gas surface density of £ = 12 
M„ pc~ 2 at t=90 Myr, the expected value from 



1 cold and dense phase with njj > 30 cm 3 and T < 300-R" 



the IKennicutt et alj (|2007l ) fit is £ S fr = 2 x 
10~ 3 M yr -1 Kpc~ 2 . This is very close to our re- 
sults. 

As observers usually do, we also calculate the 
gas consumption time-scale, r = Mqmc/SFR, in 
the simulated molecular clouds, assuming that gas 
with a density higher than 30 cm~ 3 is mainly 
within GMCs. In our simulations, the amount of 
gas in molecular clouds is Mqmc= 8 x 10 6 M at 
t=90 Myr. The star formation rate at that time is 
SFR = 4.8 x 10~ 2 M o yr" 1 . As a result, the gas 
consumption time-scale in the simulated clouds is 
r w 170 Myr. This is quite long compared with the 
typical free-fall time-scale inside molecular clouds, 
t ff = (37T/32G/9) 1 / 2 w 4 Myr for n H = 100 cm" 3 . 

In our simulations, the star formation efficiency 
over a free-fall time-scale, the fraction of gas con- 
sumed in stars during a free-fall time-scale, is only 
2.5%. This value is consist ent with observations 
( Zuckerman fc Evansl Il974h . iKrumholz fc Tanl 
(|2007j) report a range of 0.6%-2.6% for the whole 
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Kelvin (second row), gas velocity in the z-direction (third row), and surface density in cm~ 2 (forth row). 
Left panels show a face-on view of the galactic plane (z — zq) and right panels show an edge-on view 
perpendicular to that plane. The three phases of the ISM are clearly visible: cold and dense clouds, warm 
and diffuse medium and hot bubbles with very low densities. Velocities exceeding 300 km s _1 can be seen 
in hot outflows at both sides of the galactic plane. 
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population of GMCs of the Milky- Way. Our value 
is also consistent with a n efficiency of ~ 3% foun d- 



in simulations of GMCs (jVazquez-Semadeni et al 



2003tlClark et alj|2005l ). Finally, our results also 
agree with the mod el of a turbulent-dominate d 
GMC, described in iKrumholz fc McKed (120051) . 
They give an efficiency per free-fall time-scale of 
1.5%-3% for typical values of their model. 

After 100 Myr, only 10% of the gas in the sim- 
ulation has been converted into stars. Our simu- 
lations still have plenty of cold (T < 10 3 K) gas 
after 100 Myr. The surface density of this cold gas 
is ~ 5 M Q pc -2 . This value agrees with the sur- 
face density of molecular and atomic hy drogen of 



6 M pc 2 found at the solar radius (|Ferriere 



0.5 M pc 2 , compa red with the ob - 
2.5 M Q pc -2 (IFerriere I 



2001). However, the surface density of molecular 
gas is low, 
served value of 
This partially explains why our star formation ef- 
ficiency over a free-fall time-scale is in the higher 
end of the observed range. 

To conclude, stellar feedback is able to regu- 
late star formation on galactic scales because it 
regulates the amount of gas available for star for- 
mation. Stellar feedback heats and disperses the 
cold and dense gas after a star formation event. 
In a single star formation event, a stellar particle 
of ~ 90 M is created. This roughly means the 
formation of a single high-mass star embedded in 
a small stellar cluster. Due to the resolution limit, 
our simulation can not follow the details of the 
star formation process bellow ~ 10 pc scales, only 
the overall net effect. This effect is the formation 
of a small stellar cluster with an efficiency of 5%. 
As we pointed in §2.1, the star formation efficiency 
of Galactic stellar clusters is high, regardless the 
details of their formation. However, although the 
star formation efficiency is high, subsequent feed- 
back processes produce a low average star forma- 
tion. 

4.4. Volume filling factors in the ISM 

Figure [3] also shows that the net effect of stel- 
lar feedback is to produce the hot phase of the 
ISM. After the initial strong burst of star forma- 
tion, This phase can cover up to 80% of the total 
volume. This represents almost the entire volume 
above a height of 400 pc from the galactic plane. 
However, pockets of warm gas are embedded in 
this hot flow even at 2 Kpc away from the plane. 



It has the same inhomogeneous structure seen in 
the galactic fountain of figure [2] After 100 Myr, 
~ 25% of the gas is able to scape the computa- 
tional volume. 

Most of the hot gas is lost or cooled down af- 
ter 100 Myr. As a result, the volume of hot gas 
decreases because the star formation is low and 
the injection of energy is lower than in the initial 
burst. The simulation settles into a more quies- 
cent regime in which the volume occupied by the 
warm and hot phases oscillates in a self-regulated 
gas cycle. In this cycle, bursts of star formation 
(much smaller than the initial one) produce super- 
bubbles and galactic chimneys of hot gas. There- 
fore, the volume of hot gas increases. As the star 
formation fades, the bubbles cool down and the 
fraction of hot gas decreases until the next stel- 
lar burst. This pattern reflects the star formation 
history. The particular fraction of hot and warm 
phases at any moment does depend to the particu- 
lar star formation history of 10-40 Myr before that 
moment. 

4.5. Late stages of evolution 

The latter stages of the simulation offer a more 
representative view of the ISM. The effect of the 
initial conditions is gone. Therefore, we can study 
the characteristics of this feedback-driven ISM. 
We select a snapshot at 113 Myr, after the sec- 
ond burst of star formation. At that moment, 
the warm phase covers ~ 60% of the volume 
and the hot phase filled ~ 40%. X-ray emitting 
gas with temperatures above 10 5-5 K occupies 
~ 20% of the volume inside a height of 250 pc 
above the galactic plane. This is roughly consis- 
tent with ISM simulations with a 1-pc resolution 
(|de Avillez fc Breitschwerdt Il2004), Ga lactic ISM 
models and observations ( Ferrierelll998f ). 

Figure |4] shows representative slices of the box. 
The medium is very inhomogeneous at different 
scales. Large bubbles of low density coexist with 
long filamentary structures of dense clouds. Over- 
all, the medium covers more than 6 orders of 
magnitude in density and temperature. The cold 
phase forms dense and cold clouds near the galac- 
tic plane. The warm phase fills old cooled bubbles 
and low-density clouds. Finally, the hot phase 
is present in form of hot bubbles of few hun- 
dred pc wide and Kpc-scale chimneys. The gas 
in these chimneys is flowing away from the plane 
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Fig. 5. — Distribution of the gas temperature 
at 113 Myr. The distribution has three different 
peaks corresponding to three different gas phases 
of the ISM: cold, warm and hot. Two vertical 
lines show the temperature cuts used throughout 
the paper at T = 10 3 K and T = 10 4 K. 

with velocities exceeding ±300 km s _1 . These 
bubbles even break the dense plane in hot spots 
surrounded by cold and dense shells. All this 
phenomenology associated with the hot phase is 
driven by stellar feedback. As a result, one effect 
of the stellar feedback is to sustain a three-phase 
ISM. 

The distribution of temperature clearly shows 
the three main peaks of the three phases of the 
ISM (Figure [5]). The two local minima correspond 
to thermally unstable gas. The minimum around 
10 3 K, between the peaks of the cold and warm 
phases, is produced by the competition of UV 
heating and radiative cooling. This corresponds 
to the unstable regime of the classical two-phase 
model of the ISM in thermal equilibrium (Cox 
2005). The dip between 10 4 and 10 5 K results 
from the peak of the cooling curve. The gas cools 
very fast at these temperatures. As a result, it 
usually appears at the interface between hot and 
warm gas. As an exception, old bubbles at these 
temperatures are present in the simulation with 
very low densities and far away from the plane. 



0.6 
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Fig. 6. — Density distribution at 113 Myr and the 
contribution of the three phases. The dotted curve 
shows the cold phase (T < 10 3 K), the dashed 
curve shows the hot phase (T > 10 4 K), and the 
dash-dotted curve shows the warm phase (10 3 K < 
T < 10 4 K). The distribution is clearly bimodal. 
The peaks correspond to the hot and warm phases. 
The cold phase dominates the high-density tail. 

So, their cooling time is very long. This temper- 
ature distribution supports the temperature cuts 
used throughout the paper to distinguish the three 
phases: 10 3 K for the cut between cold and warm 
phases and 10 4 for the warm-hot cut. To summa- 
rize, this model of the ISM reproduces the main 
properties of the temperature distribution of the 
ISM (Cox 2005) and predicts that gas with very 
high temperatures 10 7 — 10 s K exists in the cores 
of galactic chimneys. This gas occupies only 5 % 
of the total volume and have a very small surface 
density of 4 x 1O~ 6 M pc~ 2 . 

These three phases of the ISM are also clearly 
visible in the density distribution (Figure [6]). We 
use the temperature cuts defined before to see the 
contribution of the different phases. Thus, the 
hot phase dominates the low-density range, bellow 
10~ 3 cm~ 3 . The warm phase covers intermediate 
densities, and the cold phase dominates the high 
density tail above 1 cm -3 . The density distribu- 
tion of any of the three phases cannot be described 
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Fig. 7. — Distribution of gas velocity at 113 Myr. 
The curves represents the 3 gas phases as in figure 
El Cold and warm phases have moderate veloc- 
ities, mostly bellow 100 km s _1 . The hot phase 
dominates the high velocity tail of the distribu- 
tion with velocities up to 2000 km s _1 . These are 
outflows of gas which escapes the system. 



by a single lognormal distribution, as claimed in 
Wada & Norman (2001). Instead, a combination 
of several lognorm al distributions may give a be t- 
ter approximation (jRobertson fc Kravtsovll2007l ) . 

The distribution of velocities (Figure shows 
two distinct features. The warm phase contribute 
to a strong peak around 30 km s _1 . The hot phase 
dominates the high velocity tail. It has velocities 
as higher as 2000 km s . The gas with velocities 
in this tail can easily escape the system. This gas 
forms hot outflows and galactic chimneys. 

Finally, figure [S] shows the distribution of the 
Mach number, M=u/c, where u is the gas velocity 
and c is the sound speed. The distribution shows 
that 80% of the volume has supersonic motions. 
Almost all the warm phase, half of the hot phase 
and all the cold phase are supersonic flows. The 
subsonic range is dominated by the hot phase, In 
conclusion, the ISM can hold high supersonic mo- 
tions, driven by stellar feedback. 



Fig. 8. — Distribution of the Mach number (u/c). 
The curves represents the 3 gas phases as in figure 
El 80% of the gas has supersonic motions, half of 
the hot phase has subsonic motions 



4.6. Degrading resolution 

The resolution in cosmological simulations of 
galaxy formation is much lower than the simula- 
tions of the ISM presented before. So, we can won- 
der if this picture of stellar feedback can hold if the 
resolution is degraded. Therefore, the same ISM 
models have been performed with high resolution 
(14 pc) and with low resolution (60 pc). The frac- 
tion of volume filled with each gas phase is used as 
a proxy to check the global effect of stellar feed- 
back in the ISM (Figure [9]). Left panels show that 
the hot phase covers a significant volume. 

At low resolution and without runaway stars 
(top right panel), the hot gas is almost absent 
from the simulation. Small filaments are not re- 
solved and the subsequent star formation is sup- 
pressed in these areas which can be easily broken 
by stellar feedback. As a result, star formation is 
concentrated at the center of big clumps of gas. 
Stars inject energy in high density regions, so this 
energy is radiated away without any thermody- 
namical effect in the medium. 

However, if the model of runaway stars is in- 
cluded, the hot phase is recovered at low resolu- 
tion. Stars can now migrate away from high den- 



13 



: •"' i 14 pc resolution . 

j NO RUNAWAY stars 

1 /'~~\ - 
/ \ - 
» 

I /N 

II / - 
1 \ .' 

| v..--' 
i-rvrtrf-T-t-T-T-T-T-T-v-T 


— \ / 

60 pc resolution _ 
NO RUNAWAY stars 

1 s 


:- — ■*-' ; 14 pc resolution . 

■ RUNAWAY stars 

! /''~\ 

■ 

■ !\ A 

t \ / 

\ 
I 

!— . J 

l< .'* , ;"^-4--< u. h-,— 


■ <•* ^ 60 pc resolution _ 

\ RUNAWAY stars 
\ _ 

\ 

- \r'"\/ 

V v - 
/\ / \ 
/ \ / v 
' v _ >. _--- / ' 

/ 

/ 

/ 

/ 



20 40 60 40 80 120 

t(Myr) t(Myr) 

Fig. 9. — The panels show the evolution of the 
volume occupied by each phase of the gas in four 
different models. The curves represents the 3 gas 
phases as in figure [3l The hot phase is almost lost 
for the low resolution run without runaway stars 
(top right panel). If runaway stars are included, 
the hot phase is recovered at low resolution. As 
a result, a fraction of runaway stars produces an 
effect on the global ISM and its more evident in 
low resolution runs. 

sity regions, so the injection of energy is more ef- 
ficient in forming hot gas. As a result, the model 
with runaway stars can reproduce the effect of stel- 
lar feedback even at a resolution of 60 pc. 

4.7. The expansion of a hot bubble 

As a example of the conditions of the overheat- 
ing regime discussed in section 2, we can ask now 
how a hot bubble develops in the first place. The 
top panel of figure \W\ shows the physical condi- 
tions of a single volume clement over 10 Myr. This 
volume develops a hot and dilute medium start- 
ing from a cold and dense phase. The gradients 
are computed using a 3-points finite differences ex- 
pression using the adjacent cells. 

At the beginning, there are no stars present in- 
side that volume. So, there is no feedback heating. 
At the same time, the density is high enough so the 
the radiative cooling dominates over the UV back- 
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Fig. 10. — Evolution of the properties of a single 
volume element (Temperature, density, gradients 
of pressure and gravity, and the mass in young 
stars) in a run with 8 pc resolution (top panel) 
and with a resolution of 30 pc (bottom panel). 
The top panel shows an over-pressured volume 
that expands due to stellar feedback. It produces a 
hot cavity filled with low-density gas. The bottom 
panel shows how this hot bubble cannot develop 
when the gradients of pressure do not dominate 
over gravity. 

ground heating. As a result, the medium stays at 
the floor temperature of 300 K. The medium is 
also in hydrostatic equilibrium. 

The situation drastically changes when young 
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stars appear. They are not born inside that partic- 
ular volume. Instead, they are drifting slowly from 
adjacent cells. The result is that this young popu- 
lation injects energy into the medium. So, heating 
dominates over cooling initially. The system re- 
sponds by increasing the temperature. As a result, 
the cooling rate increases and the medium reaches 
a balance between cooling and heating rates in a 
very short time-scale. This is because the cooling 
time is very short in those conditions. The net 
result is a medium slightly hotter than surround- 
ings so this over-pressured region expands and the 
density inside that volume decreases. 

Around 10 4 K, the cooling curve is a very steep 
function of temperature, so the temperature in- 
creases very slowly. But, at the same time, the 
density drops faster. As a result, the cooling rate 
decreases. This expansion is fueled by a roughly 
constant injection of energy from massive stars. 

When the conditions of eq. are fulfilled, the 
medium can pass through the peak of the cooling 
curve, somewhere between 10 4 — 10 5 K. After that, 
the gas has low density and a temperature of few 
millions degrees. As a result, heating dominates 
over cooling and a hot cavity is formed. 

The bottom panel of figure [TOl shows a different 
situation. The volume is selected to be the high- 
density core of a molecular cloud formed in the 
low-resolution run shown in the top-right panel of 
figure [5] Stellar feedback from the stars formed in 
that core are able to heat the gas only to 10 4 K. 
The gradients of pressure do not overcome gravity. 
The condition of bubble expansion is not fulfilled , 
eq. (HI). As a result, the density remains high and 
a hot bubble can not develop. 

5. Results of Cosmological runs 

In previous sections, we have shown that our 
models of stellar feedback follow the effect of su- 
pernova explosions and stellar winds in the ISM 
with a resolution of about 50 pc. The result is the 
formation of super-bubbles and galactic chimneys. 
Both arc filled with hot and dilute gas. The net 
result is a multi-phase ISM and galactic outflows 
with large velocities. 

Now, we can study the effect of stellar feed- 
back in galaxy formation. We apply these feed- 
back models in cosmological hydrodynamics sim- 
ulations with a similar resolution of 35-70 pc. The 




Fig. 11. — Density- weighted average gas density 
along the line-of-sight through a MW progenitor 
at redshift 3. The top panel shows a smooth gas 
distribution with some density enhancement due 
to a pattern of spiral arms. Young stars appear 
as points. They follow the spiral pattern. How- 
ever, the distribution has a very concentrated and 
dense center. The bottom panel shows the heating 
regime case. The distribution is not as dense, so 
the range of density is different. The center is less 
concentrated and the distribution looks clumpy 
with dense filaments and clouds embedded in a 
more diffuse medium. This is the case of a multi- 
phase ISM. The size of the images is 30 Kpc in 
both cases. 

simulations follows the formation of a MW-type 
galaxy starting from primordial density fluctua- 
tions. 

The computational box is 10 hr 1 Mpc com- 
moving box. We apply a zooming technique 
( Klypin et al. 200lh to select a lagrangian vol- 
ume of 3 virial radius centered in a MW-size halo 
at redshift 0. Then, we resimulate that volume 
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Fig. 12. — Circular velocity profile of the main 
progenitor of a MW-type galaxy at redshift 5. The 
top panel shows the results of an inefficient stellar 
feedback. The galaxy is too concentrated and has 
a too massive spheroidal component. By contrast, 
the bottom panel shows a regime in which stellar 
feedback is more efficient and it can regulate the 
growth of the galaxy. The maximum resolution in 
both cases is 70 pc. 



with higher resolution. The region has a radius 
of about 1.5 hT 1 comoving Mpc. The simulation 
has about 5 million dark matter particles. They 
have three different masses. The high-resolution 
region is resolved with 3.4 million dark matter 
particles with a 7.5 x 10 5 M mass per particle. 
The high-resolution volume is resolved with about 
17 million volume elements at different levels of 
resolution. The maximum resolution is always 
between 35 and 70 pc. A short summary of the 
details of the simulations is given in table [TJ The 
cosmological model assumed throughout the paper 
has tt m = 0.3, Qa = 0.7, h=0.7 and a s = 0.9. 

5.1. Heating regime versus overcooling 
regime 

We compare two cosmological simulations with 
the same spatial resolution but different regimes. 
Table [1] shows a summary of the two simulations. 
The over-cooling model has low star formation ef- 



Fig. 13. — Distribution of velocities at the highest 
levels of the resolution for two models. The dashed 
curve represents the hot phase with T > 10 4 K and 
the dash-dotted curve represents gas with temper- 
atures bellow 10 4 K. The bottom panel shows a 
longer tail of high-velocities. These are galactic 
outflows which can reach velocities exceeding 10 3 
km s . 

ficiency. In addition the cosmological UV ba ck- 
ground according to lHaardt fc Madau ] (|l996h is 
present over the whole evolution. Finally, a con- 
stant model of stellar feedback is used. 

In the second simulation the heating regime de- 
velops. It has a high efficiency of star formation 
and the UV background is limited to its value at 
redshift 8. The supernova model of stellar feed- 
back was used in this case. We use a model of star 
formation in which each star formation event was 
treated as a random process (see Appendix). In 
this way, we keep a moderate galactic star forma- 
tion rate of ~ 10 M yr _1 inside the main galaxy 
at redshift 3. 

The simulation in the cooling regime has a cold 
galactic ISM with temperature close to 10 4 K. 
The simulation in the heating regime develops a 3- 
phase ISM. Hot bubbles develop naturally. They 
produce galactic chimneys that combine in a galac- 
tic wind. As a result, galactic winds are the natu- 
ral outcome from stellar feedback. 
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Table 1 

Parameters of cosmological models. 



Parameter 



Models 



Comoving box size 
Number of DM particles 
DM particle mass 
Number of cells 
Max. resolution (proper) 
Max. number of stars 
Min. mass of stellar particle 
Model name 
UV flux 

Star formation time scale r 
Model for stellar energy release 
Runaway stars 



Overcooling 
H& M96 

4 x 10 7 yrs 
Constant 
not included 



14.28 Mpc 
5.4 x 10 6 
7.5 x 1O 5 M 
17.5 x 10 6 
35-70 pc 
3.7 x 10 6 
10 4 M Q 

Multi-phase 
H& M96 but constant after z 
4 x 10 6 yrs 
SN model 
included 



Figure [TT] shows the main MW progenitor at 
redshift 3 for both simulations. The cooling model 
has a smooth density distribution with a small en- 
hancement due to a pattern of spiral arms. In 
contrast, the multi-phase model develops a clumpy 
medium of dense clouds surrounded by low-density 
bubbles. This is a multi-phase medium. 

5.2. Comparison of circular velocity pro- 
flies 

We can see now the effect of this multiphase 
medium in the galaxy assembly. We use the pro- 
file of circular velocity as a proxy of the mass dis- 
tribution, V c — y/GM/R, where G is the gravita- 
tional constant and M is the mass inside a radius 
R. Figure [12] shows the profile of the circular ve- 
locity for the same galaxy in the two cases. The 
simulation with the overcooling problem shows a 
strong peak in the baryonic component of the cir- 
cular velocity. Both gas and stars are very concen- 
trated in the first Kpc. In contrast, the simulation 
with multiphase medium has a more shallow cir- 
cular velocity profile. This indicates a less con- 
centrated galaxy with less baryons. At the virial 
radius, K V i r —16 Kpc, the virial mass at z—5 is 3.1 
xlO 10 M . A large fraction of this mass is dark 
matter, M dm = 2.7 xlO 10 M . The mass in gas 
is M gas = 0.24 xlO 10 M Q and the baryons locked 
into stars accounts for only 0.14 xlO 10 M . 



5.3. Galactic winds and multiphase medium 

The hot bubbles in the multi-phase medium de- 
velop galactic fountains that produce hot outflows 
with very high velocities: larger than 10 3 km s _1 . 
These outflows are not produced in the cooling 
model. The Figure 1X31 shows the difference in the 
distribution function of velocities for both cases. 
We take all cells at the highest levels of refinement. 
Therefore, we select a volume close to the galaxies 
in the simulations. The multi-phase model has a 
bigger fraction of hot gas with much larger veloc- 
ities than in the cooling model. These outflows 
contribute to the high- velocity tail of the distribu- 
tion. In the cooling model, the distribution drops 
at 300 km s _1 , while in the hot case the tail ex- 
tends beyond 10 3 km s _1 . 

These galactic-scale outflows can be seen figure 
[Ml It shows a slice of the simulation through the 
main MW-progenitor at redshift z = 3.4. At that 
redshift, its virial radius is 70 Kpc and the total 
virial mass is 10 11 M©. The gas density panel 
shows the galaxy embedded in a cosmological web 
of filaments. The galaxy at the center is blowing 
a galactic wind of hot and dilute gas with out- 
flows velocities exceeding 300 km s _1 . The wind is 
rich in a-elements and other products of the ejecta 
of core-collapse supernova. These metal-rich out- 
flows can contribute to the enrichment of the halo 
and the inter-galactic medium. These outflows can 
reach even higher velocities and can escape the 
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Fig. 14. — This panel shows a galaxy at redshift 3.4 with a resolution of 45 pc. The figures show slices 
of 600 Kpc on a side of gas density (top left), temperature (top right), velocity in the horizontal direction 
(bottom left), and metallicity (bottom right). There are inflows of low-metallicity gas in cold filaments, as 
well as outflows of hot, metal-rich gas produced by chimneys in a multi-phase interstellar medium. Outflow 
velocities exceeds 300 km s _1 . The virial radius is 70 Kpc and the total virial mass is 10 11 solar masses. 




Fig. 15. — The same as in figure [T4l but now the size of the images is 50 Kpc. It shows a multi-phase ISM 
of cold and dense clouds surrounded by bubbles of hot and dilute gas. Inflow and outflows velocities can 
reach 300 km s _1 . The outflows are galactic chimneys powered by core-collapse supernova. Therefore, they 
are rich in a-elements. In contrast, the inflow of gas has almost primordial composition. 
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R (Kpc) 

Fig. 16. — Density profile of a galactic halo at red- 
shift 5. The dark matter distribution is consistent 
with a cusp. The diagonal line shows a slope of -1. 

galactic halo and enrich the inter-galactic medium. 
The galactic wind is produced by the combina- 
tion of different galactic chimneys anchored in the 
multi-phase ISM of the galaxy. 

Figure [15] shows this multi-phase ISM. Cold 
and dense clouds coexist with low-density bubbles 
filled with very hot gas. Warm gas with inter- 
mediate densities and temperatures filled areas of 
low star formation and inflows of gas with almost 
primordial composition. 

5.4. Density profile consistent with a 
cuspy profile 

Figure [16] shows the inner profile of density 
of the different components of the galaxy: dark 
matter, gas and stars at redshift 5. The density 
slope of the dark matter profi le is consistent with 
a cusp y profile. In contrast, iMashchenko et a"l~l 
(|2007l ) reported the formation of a core rather than 
a cusp in the central ~ 300 pc of a much smaller 
galaxy at high redshift ( ~ 10 9 M at z=6) in a 
SPH cosmological simulation. In their case, the 
mechanism that removes the cusp is gravitational 
heating from large fluctuations in the gravitational 
field. These fluctuations are produced by bulk 
motions of gas clumps driven by stellar feedback 



(jMashchenko et al.|[2006h . These motions remove 
episodically 90% of the mass from the central 100 
pc after each burst of star formation. 

However, these gas clumps can be overproduced 
in simula tions if the local Jea ns length is not 
resolved (jTruelove et al. 1 119971 ). This produces 



an artificial gas fragmentation and big clumps of 
stars. An excessive dumpiness can artificially in- 
crease the efficiency of this gravitational heating. 
In our simulations, we prevent this artificial frag- 
mentation by the implementation of a pressure 
floor that increases the effective Jeans length t o 
the resolution limit ( Robertson fc Kravtsovi2007l ). 



Ho wever, a direct comp a rison between our results 
and Mashchenko et al. (2007) is difficult because 
we follow the formation of a much bigger galaxy (~ 
10 10 M Q at z=5), in which the effect of this grav- 
itational heating driven by stellar feedback is less 
important. Therefore, these gravitational heating 
driven by stellar feedback can not be ruled out in 
low-mass and gas-rich starburst galaxies. 

6. Summary and conclusions 

We study the role of supernova explosions and 
stellar winds in the formation of galaxies. Our ap- 
proach is to model these processes without the ad- 
hoc assumptions typically used on stellar feedback. 
Unlike many currently used prescriptions, we do 
not stop cooling in regions where the energy from 
stellar feedback is released (Thacker & Couchman 
2000; Brook et al. 2004; Keres et al. 2005; Gov- 
ernato et al. 2007). Moreover, instead of using a 
sub-resolution model of a multi-phase medium ( 
Springel & Hernquist 2003, Cox et al. 2006), we 
resolve that multi-phase medium. This is a more 
straightforward way to model stellar feedback. It 
eliminates many ad-hoc assumptions. This ap- 
proach also produces naturally the outcomes usu- 
ally associated to stellar feedback: hot bubbles, 
chimneys and galactic winds. 

Feedback heating has an effect in the ISM only 
when it dominates over radiative cooling. Section 
2 shows the necessary conditions for this heating 
regime (eq.(|5|[7])). We find that a model of cooling 
bellow 10 4 K is a key ingredient to fulfill these 
conditions. Thus, by resolving the conditions of 
molecular clouds (T « 100 K and nn > 10 cm -3 ), 
we resolve the conditions, in which stellar feedback 
is more efficient in the ISM on galactic scales. 
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We perform parsec-resolution simulations of a 
piece of a galactic disk in order to see the effects 
of stellar feedback and to test our models. When 
we use a realistic feedback and high resolution, 
the system has a low star formation rate and it 
forms hot super-bubbles of 100-pc scales and Kpc- 
scale galactic chimneys. We found that the cores 
of these chimneys reach temperatures of 10 7 — 10 s 
K, very low densities [hh < 10~ 4 cm -3 ), and out- 
flow velocities exceeding 10 3 km s -1 . 

Then, we degrade the resolution to see if this 
picture of multi-phase ISM holds at a resolution 
that we can achieve in cosmological simulations. 
We found that runaway stars help to spread the 
effect of stellar feedback. They usually explode as 
supernovae in low-density regions, few 100 pc away 
from their natal molecular cloud. This is an effect 
found in nature (Stone 1991), which enhances the 
feedback. So, it should be included in any realistic 
model of stellar feedback. 

Thermal feedback from young stars is able to 
produce long time-scales of gas consumption by 
dissipating the star-forming gas. As a result, 
although this gas has high star formation effi- 
ciency, subsequent feedback processes produce a 
low star formation rate, averaged over all cold 
and dense gas. For example, in the simulations 
of the ISM described in §4, the gas with a den- 
sity above the density threshold for star forma- 
tion can form stars with high efficiency. How- 
ever, the average star formation efficiency in the 
simulated clouds is roughly 2.5% over a free- 
fall time-scale (§4.3). This is roughly consis- 
tent with estimations of the star formation effi- 



ciency in molecular cloud s (Zuckerman & Evans 



1974; lKrumholz fc McK ce 2005 UKrumholz fe Tan 



20071) . 

In cosmological simulations (§5), We find a 
moderate galaxy star formation rate, SFR = 10 
M yr _1 and a significant amount of cold and 
dense star-forming gas, Mdensegas = 1O 9 M in- 
side a 5 Kpc star- forming disk at redshift 3. These 
values are consistent with observations of nearby 
starburst galaxies. Using the observed relation 
between the star formation rate and the amount 



of star- forming gas of lGao fc Solomon! (|2004f ). the 



star formation rate expected for 10 M of cold 
and dense gas is 20 M yr _1 . This is close to the 
value found in your simulations. Moreover, the 
galactic gas consumption time-scale of dense gas, 



Mdcnsogas/SFR is ~ 100 Myr. This is consistent 
with observed values in local starburst galaxies 
which are usually used a s analogs of star-f orming 
galaxies at high redshift ( Kennicutt l fl998h . 

In our simulations, star formation proceeds in 
a way co nsistent with obser vations of star- forming 



galaxies ( Kennicutt 1998). From the numbers 



given above, the gas surface density of the star- 
forming disk of 5 Kpc radius at redshift 3 is 
S gas = 13 M pc~ 2 . Using th e Kennicutt fit fo r 
nearby star-forming galaxies (jKennicutt 1998). 



the expected value for the star formation rate sur- 
face density is E SFR = 10~ 2 M yr -1 Kpc" 2 . The 
measured value from the simulations is Ssfr = 
1.3 x 10 _1 M yr^ 1 Kpc -2 . Although this value is 
an order of magnitude higher than the expected 
value from the fit, it is still within the intrinsic 
spread found in observations. As a result, our 
simulated high-redshift galaxy seems more com- 
pact than the average star-forming galaxy at low- 
redshift. 

Our cosmological simulations with this model of 
stellar feedback do not have the overcooling prob- 
lem. The fraction of cold baryons (stars and gas 
with a temperature bellow 10 4 K) inside the virial 
radius at z—5 is 0.6 times the cosmological value ( 
■fr.o.smn= 0-15). This is cons istent with galaxy mass 
models (jKlvpin et al 1l2002h . Instead of a cold disk, 
we produce a multi-phase ISM with the same fea- 
tures seen in the simulations of the ISM described 
in section 4: cold clouds, hot super-bubbles and 
galactic chimneys. The angular momentum prob- 
lem is also reduced. Instead of a compact object 
with a strong peak in the rotation curve, we pro- 
duce more extended galaxies with nearly flat rota- 
tion curves. Baryons are less concentrated when 
stellar feedback plays a role in the formation of 
galaxies. At the same time, the density profile of 
dark matter is still consistent with a cuspy profile. 

In this picture, galactic chimneys powered by 
stellar feedback combine into a galactic wind. So, 
galactic winds appear as the natural outcome of 
stellar feedback in starburst galaxies at high red- 
shifts. We found typical outflow velocities of 300 
km s _1 with some exceptional examples of out- 
flows exceeding 1000-2000 km s _1 . This is consis- 
tent with observa tion of outflows at high redshift 
(|Law et al. II2007I ) . From a s ample of ~ 100 galax- 



ies at redshift 1.9 < z < 2.6. ISteidel et al. 



find a mean outflow velocity of 445 km s 



1 2007ft 

Some 
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cases have velocities of 1000 km s . 

This picture is only reproduced if the resolu- 
tion is high enough to resolve the physical condi- 
tions of densities and temperatures of molecular 
clouds. Our cosmological simulations reach a res- 
olution of 35 pc, which is 10 times better than the 
typical resolution in previous cosmological simu- 
lations (Sommer-Larsen et al. 2003; Abadi et al. 
2003; Robertson et al. 2004; Brook et al. 2004; 
Okamoto et al. 2005; Governato et al. 2007). 

Acknowledgments 

We are grateful to A. Kravtsov for providing the 
hydro code. We are in debt to N. Gnedin for giv- 
ing us the wonderful analysis and graphics package 
IFRIT. We thank K. Tassis for very useful dis- 
cussions. We acknowledge support of NSF grants 
to NMSU. The computer simulations presented in 
this paper were performed at the National Energy 
Research Scientific Computing Center (NERSC) 
of the Lawrence Berkeley National Laboratory and 
at the NASA Advanced Supercomputing (NAS) 
Division of NASA Ames Research Center. 



21 



A. A model of star formation for scales bellow 100 pc 

A successful model of star formation in simulations should take into account the spatial resolution. For 
example, in typical cosmological simulations with a resolution of ~ 1 Kpc, the star formation is averaged over 
a large piece of ISM. These simulations should have a star formation model with a low star formation efficiency 
in order to reproduce the global efficiencies found in nearby galaxies. Observations of quiescent galactic disks 
show long gas consumption time-scales averaged over a significant piece of a galaxy, T" g i ODa i = Ep-a. s/^RFR ~ 1 



Gyr , where Esfr is the st ar formation rate surface density an £ gas is the gas surface density ([Kennicutt 



ll998tlKennicutt et al.1120071 ). At t he same time, for s tarburst galaxies, the global gas consumption time-scale 



is much shorter, r g i b a i = 0.1 Gyr ([Kennicutt 1119981 ) 



However, if the resolution is high enough to resolve the regions where star formation mainly occurs, 
giant molecular clouds, the star formation efficiency can be much higher: the time-scales for the formation 



of Galactic stellar clusters a re around few Myr and 10% - 40% of the gas is consumed (jGreene fe Young 



1992; Elmegreen et al.ll2000h . As a result, simulations which can resolve the sites of star formation should 



have a high star form ation efficiency only in the high-density regions, where molecular clouds can form 
(|Tasker fc Brvanll2006l) . [n practice, the maximum resolution that we can afford is between 30-70 pc. This 



limits the maximum densit y that our simulations c an resolve. For example, if we consider a typical giant 
molecular cloud of 10 5 M ( Rosolowskv et al~ll2007l ) . the mean density averaged over 30-80 pc scales will be 
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10-200 cm . This gives an idea of the typical densities where star formation occurs our simulations. 

In our code, star formation is allowed in a time step, dtsF> which is equal to the time step of the 0- 
Level of resolution. This time step is controlled by the Courant condition for hydrodynamics and in our 
cosmological simulations, dt$F — 1-2 Myr. During this period of time, a stellar particle can form only 
where the density and temperature reach a given threshold: p gas > psf and T gas < Tsf- Even in these cold 
and dense regions, each star formation event is treated as a random event with a probability Pr to occur. 
We roughly approximate the fact that regions with higher densities have a higher probability to host star 
formation events by assuming a simplified formula: 

Pr = -^5_ (Al) 
100p SF 1 ; 

In this way, the number of stellar particles remains in a value that is not computational prohibit ed. In the 



form ation of a single stellar particle, the star formation rate is proportional to the gas density (jKravtsov 
20031) : 

dp*, young _ Pgas (A2) 
dt T 

where p*, y0 ung is the density of new stars, p gas is the gas density and t is a constant star formation timescale. 
The density and temperature thresholds used are psf = 0.035 M Q pc~ 3 = 1 cm -3 ) and T$f = 10 4 K. 
In spite of the fact that we allow star formation starting at 10 4 K, in practice the vast majority (> 90%) of 
"stars" form at temperatures below 1000 K and more than half of the stars form bellow 300 K and densities 
larger than 10 cm~ 3 . 

As described in §2.1, the ratio p*, y0 ung/Pga,s should be ~ 0.1-0.5 for typical conditions of dense, star- 
forming gas. Only in this case thermal feedback can produce over-pressured hot bubbles in the sites of star 
formation (eq. [5]). Based on equation IA21 this ratio of densities can be expressed as 

P*, young _ ^SF (A3) 



As a result, thermal feedback is only efficient in dense, cold, star-forming gas if dtsF/T ~ 0.1 — 0.5. This sets 
the value of r, because dtsF is set by the conditions of hydrodynamics, as explained before: dtgF = 1 — 2 Myr. 
Therefore, the value of r should be in the range 2-20 Myr, consistent with the gas consumption time-scales 
during the formation of Galactic stellar clusters ( Greene fc Younel 1992 : Elmegreen et al. I l2000h . However, 
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this high local efficiency of star formation in high-density regions produces the observed low global efficiency 
Tgiobai = 0.1 — 1 Gyr, as discussed in §4.3. 
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